######################
##### PREPARATION ####
######################
# Clear environment
rm(list=setdiff(ls(), "path_to"))

# Load data
s = readRDS(path_to("wide"))


######################
#### DERIVE DATA #####
######################

# Distribution
factor_levels = levels(s$wtp_ratio_class)
s$wtp_ratio_class = as.character(s$wtp_ratio_class)
s$wtp_ratio_class[is.na(s$wtp_ratio_class)] = "not defined"
dist = rbind(
  s %>% group_by(conseq, wtp_ratio_class) %>% summarise(
      n = n(),
      wgt = sum(wgt),
      nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
      consumer = sum(consumer_vig),
      dampening = sum(belief_dampening == 1)
    ) %>%
    mutate(
      freq_main = n / sum(n, na.rm=T),
      freq_wgt = wgt / sum(wgt, na.rm=T),
      freq_nospeeders = nospeeder / sum(nospeeder),
      freq_consumer = consumer / sum(consumer),
      freq_dampening = dampening / sum(dampening),
      wtp_ratio_class=factor(wtp_ratio_class)
  ),
  s %>% group_by(wtp_ratio_class) %>%
    summarise(
      n = n(),
      wgt = sum(wgt),
      nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
      consumer = sum(consumer_vig),
      dampening = sum(belief_dampening == 1)
    ) %>%
    mutate(
      conseq = "pooled",
      freq_main = n / sum(n, na.rm=T),
      freq_wgt = wgt / sum(wgt, na.rm=T),
      freq_nospeeders = nospeeder / sum(nospeeder),
      freq_consumer = consumer / sum(consumer),
      freq_dampening = dampening / sum(dampening),
      wtp_ratio_class=factor(wtp_ratio_class)
    )
)
dist = dist %>% pivot_longer(starts_with("freq_"), values_to="freq")
dist$wtp_ratio_class = factor(dist$wtp_ratio_class, levels=c("not defined", factor_levels))
dist$name = gsub("^freq_", "", dist$name)

# Order
conseqs = c(
  "pooled" = "All cases",
  "co2" = "CO2 emissions\n(Reduce by 1 ton)",
  "waste" = "Non-recyclable waste\n(Reduce by 100 pounds)",
  "chicken" = "Animal welfare\n(20 certified chicken)",
  "textile" = "Low wages\n(10 fairtrade garments)"
)
dist$conseq = factor(dist$conseq, levels=names(conseqs))

# Robustness tests
tests = c(
  "main" = 1,
  "wgt" = 2,
  "nospeeders" = 3,
  "consumer" = 4,
  "dampening" = 5
)


######################
#### PLOT ############
######################
plots = subplots = list()

plots$fig_conseq_robust = ggplot(dist, aes(x = name, y = freq, fill = wtp_ratio_class)) +
  geom_bar(stat = "identity") +
  facet_grid(.~conseq, labeller = as_labeller(conseqs)) +
  xlab("Specification") + ylab("Frequency") +
  scale_x_discrete(limits = names(tests), labels=tests) +
  scale_y_continuous(breaks=seq(0, 1, 0.2), labels=paste0(100*seq(0, 1, 0.2), "%")) +
  scale_fill_manual(
    values=c("gray90", "red3", "#e06258", "#d68f8e", "#cfb5ba", "#b9c2d0", "lightblue2"),
    labels=paste0(c("(Not valued)", factor_levels), " "),
    name="**Valuation ratio**: valuation if ineffective / valuation if effective *(from left to right)*"
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "rob_con"),
    plot.title = element_text(size=18, face="bold", hjust=0.5),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.title = element_markdown(size=14, hjust=0.5),
    #legend.spacing.y = unit(0, 'cm'),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    strip.text = element_text(size=14, face="bold"),
    axis.title = element_text(size=14, face="bold"),
    axis.text.x = element_markdown(size=14),
    axis.text.y = element_text(size=14),
    legend.text = element_text(size=14)
  ) +
    guides(fill = guide_legend(nrow = 1))

# Save generated graph.
pdf(path_to("figures", "fig_conseq_robust.pdf"), width=11, height=5.5)
print(plots$fig_conseq_robust)
dev.off()  


